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1. Introduction 

The airburst on February 15, 2013 of a superbolide (approximately 20 m in size) over 
the Russian city of Chelyabinsk has given a new impetus, and a renewed sense of 
urgency, to planetary defense programs across the globe. In addition to the ongoing 
efforts in the detection and characterization of Potentially Hazardous Asteroids 
(PHAs), and developing strategies to mitigate the risk of potentially devastating 
impacts of large size (0.5 km or larger) PHAs, there is new interest in understanding 
and quantifying the risk of airburst/impact of PHAs that are much smaller than 0.5 
km. 

In October 2014 NASA’s Near Earth Object (NEO) Program Office initiated a new 
effort in Planetary Defense (PD) at Ames Research Center (ARC). A companion to 
the present paper, by Arnold and Burkhard [1], provides the objectives and details of 
this new program, which integrates experiment and analysis in the areas of planetary 
science [2], entry physics (present work), and blast wave physics into a physics- 
based risk assessment framework [3]. Briefly, the research endeavor aims to 
determine scenario-based probabilities of the outcomes of impacts of PHAs of 
various sizes and various spectral classes (stony, stony-iron, iron, etc.), thus helping 
decision makers develop mitigation strategies. The present paper provides an 
overview of the work conducted for the entry physics task of the ARC PD program. 


2. Objectives 


We have undertaken the effort to extend and apply tools used for design and 
analysis of atmospheric entry capsules to the problem of entry of meteoroids, with 
the overarching objectives of understanding energy deposition into the atmosphere 
and fragmentation . Since the effort is still in its early stage, the near term objectives 
(for the purposes of the present paper) are; (i) to explore the physics of atmospheric 
entries of meteoroids using our current state-of-the-art tools and processes, (ii) to 
verify results of computations against results from models widely accepted by the 
meteor physics community, (iii) to explore the influence of shape (and shape 
change) on flow characteristics, and (iv) to explore how multiple bodies interact. In 
time we hope to improve both the tools and our understanding of the physics of 
energy deposition and break up, thereby adding more fidelity into the risk 
assessment tools. 

3. Methodology and Assumptions 

The methodology we have chosen to adopt in the initial phase of our effort is briefly 
described below. 

1 . We parameterize the flight space (the concept is described in more detail in the 
work of Prabhu et al. [4]) by the flight velocity, freestream density (altitude), and 
object size in an attempt to determine scaling laws for use in trajectory 
computations. This parameterization delinks flight trajectories from computations 
of the aerothermal environments around the meteoroid. 

2. We consider only regular (and smooth) geometries, such as spheres and prolate 
spheroids. Having one axis (or axes) of symmetry considerably simplifies the 
effort in generating volume meshes. Furthermore, the use of simple geometries 
helps us tie into the work done in the meteor physics community. We note that 
the flow computations are not necessarily restricted to single objects. We also 
perform 3D computations for a finite collection of objects (still regular smooth 
geometries) to understand the interactions of the objects. 

3. Flow computations are performed in the parameterized flight space without 
making assumptions about the meteoroid material and/or internal structure. 
Computations include estimation of high-temperature shock layer radiation. The 
results of computations at the end of this step are to be considered bounding 
because they do not include the thermal response of the material being heated 
by convection and radiation. 

4. The computed environments are then used in: (1) trajectory software 
development, (2) material structural and thermal-structural response, and (3) 
material thermal response. The expectation is that the responses will be fed back 
into flow simulations, i.e., the thermal and structural response models are loosely 
coupled to the flow model. The feedback (chemical species and shape change) 
from the response models helps quantify the amount of energy attenuation due to 
ablation. 

The altitude-velocity plot shown in Fig. 1 contains the flight space (quadrilateral) 
used in the present work. Also shown in the figure are the best-estimated trajectory 
(BET) for Stardust [5], and data for the Chelyabinsk bolide [6]. At each point shown 
in the flight space (filled diamonds in Fig. 1), flow computations are performed for 
regular geometries (hemispheres for now) of sizes ranging from 1 m to 100 m. 
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Figure 1. Flight space parameterization for flow computations. A velocity range of 12 to 20 
km/s is covered for stagnation pressures ranging from 0.1 to 300 bar (0.01-30 MPa). 

As is seen from Fig. 1, we have considered velocities from 12 to 20 km/s 
(equivalently 72 MJ/kg to 200 MJ/kg of kinetic energy). The state-of-the-art analysis 
tools and processes are calibrated to the Stardust mission (see the special issue of 
Journal of Spacecraft and Rockets, 47(6), 2010), which is at lowest end of the 
velocity space that we have to explore for meteor entries. Therefore, application of 
these tools to Chelyabinsk-like meteor entries is a significant challenge. 

From the computed flowfields, the surface pressure, shear stress (magnitude), and 
heat flux are extracted. The drag coefficient, Cd, is determined by integrating the 
pressure and shear stress distributions, and the heat transfer coefficient, Ch, is 
determined by integration of the total heat flux (convective and radiative) 
distributions. Finally, with a view towards developing a synthetic light curve, the 
energy radiated outward (/.e., luminosity) from the shock layer into a hemisphere 
enclosing the bow shock is calculated. Although radiation energy flux computations 
are performed over the wavelength range of 85-39600 nm (VUV to IR), the 
wavelength range is restricted to 200-800 nm (middle UV to visible) in the 
computation of luminosity. 

4. Analysis Tools 

The primary modeling and simulation tools are: 

1. Traj [7]: A 3DoF trajectory simulation code. The code has a variety of 
atmospheric models (for all relevant planets and moons in the Solar System), and 
gravitational and rotating planet models. For the present work, it has been 
updated to include a simple differential equation for mass loss, with heat transfer 
coefficients obtained from detailed flowfield computations via correlations. 

2. Dplr [8]: A 3D Navier-Stokes flow solver for hypersonic flows. The code includes 
models for shock-layer thermochemical nonequilibrium, and a variety of surface 
boundary conditions. In the present work we have an 11-species air model (N 2 , 


02, NO, N 2 '', O 2 '', NO"", N, O, N"", O'", and e') in thermal equilibrium (i.e., a single 
temperature is used for the flowfield). This choice of model limits us to velocities 
not exceeding 20 km/s, since thermodynamic, transport, and rate data for higher 
stages of ionization of N and O are still being developed. In all computations the 
wall (surface of the meteoroid) is assumed to be at 300 K and fully catalytic to 
atom recombination. Radiative equilibrium (hot wall) computations have been 
deferred until the material thermal response code is fully operational for 
meteoritic materials. 

3. Neqair [9]: A line-by-line spectral code with the tangent-slab approximation for 
transport of radiation energy. The code is applied to lines of sight distributed over 
the surface. The data along a line of sight, T and ns (temperature and number 
densities of species) are extracted from the flow field solution. The methodology 
developed by Palmer et al. [10] for loosely coupling the radiation code, Neqair, to 
the flow solver, Dplr, has been applied to a small number of cases in the flight 
space. Since such a coupled computation is not only expensive but also 
somewhat sensitive to the frequency of updates/interchanges between the flow 
and radiation fields, an empirical relation [11] is used to correct the adiabatic 
value obtained by a posteriori computations of radiation from the other converged 
flow field solutions. 

4. Fiat [12] and Titan [13]: One- and two-dimensional material thermal response 
codes. These codes rely on the idea of vapor blowing from the heated surface. 
The non-dimensional blowing rates, referred to as B-prime tables, are 
constructed separately for the desired range of surface temperatures and 
pressures. B-prime tables are readily available for materials typically used in the 
thermal protection system of entry capsules. However, construction of such 
tables for each meteor type (spectral class) is ongoing work. These codes too 
can be loosely coupled to the flow solver, Dplr. 

5. Marc [14]: Commercial finite-element analysis (FEA) software for structural and 
thermal-structural response computations. The software package supports fully 
transient, nonlinear, coupled thermal-structural FEA, and can use nonlinear and 
temperature-dependent material properties. The package allows the application 
of a wide range of structural and thermal boundary conditions that include applied 
pressure, shear, displacement and rotation constraints, heat flux, convective 
heating, fixed temperature, and radiation. The boundary conditions can be simple 
constant values, or they can vary both spatially and temporally. 

6. Utility programs [15]: These are numerous CFD-related programs developed in- 
house at NASA ARC for various projects (past and ongoing). 

5. Results and Discussion 

We next present a sampling of results from what we have accomplished thus far. 
The results are from exploratory studies and are best considered preliminary, 
especially since models for: (a) response of meteoritic material to imposed heating, 
and (b) structural and thermal-structural analysis are still under development while 
relevant properties (at higher temperatures) are being gathered. 

5.1. Trajectory simulations 

The mass loss equation used in meteor physics has been added into the trajectory 
simulation tool, Traj [7]. The code currently allows for the specification of a heat 


transfer coefficient as a constant over the entire trajectory. Modeling of dependence 
of this coefficient on the velocity, altitude (density) and size of a spherical geometry 
is still under development. However, in the interim, the trajectory tool has been used 
to determine the shallowest possible entry flight path angle over ballistic coefficients 
ranging from 0 to about 500,000 kg/m^ (Fig. 2). 



Skip out angle/deg 

Figure 2. Skip out angle for ballistic coeffcients ranging from 0 to 500,000 kg/m^ and four entry 
velocities. 

These results are useful for exploring the influence of entry flight path angle (defined 
as the angle reckoned from the local horizon with the negative sign indicating that 
the velocity is below the local horizon). A simple definition of the entry ballistic 
coefficient (for a spherical geometry) is yS = {4I3)pR/Cd, where p is the bulk density of 
the meteoroid, R is the spherical radius, and Cd is the hypersonic drag coefficient of 
a sphere (»0.9). Therefore, for a given material class (bulk density) and size, the 
limiting entry angle on the shallow end can be simply read off Fig. 2. 

5.2. Aerothermal environments database for hemispheres 

Axisymmetric flow computations were performed on hemispherical geometries (i.e., 
a wake is not included) for diameters ranging from 1 m to 100 m and for velocities 
ranging from 12 to 20 km/s. The objective was to determine the heating 
environments (convective and radiative) and visible light emission from the gas cap. 
In all cases the material is passive (no ablation) and the flow is assumed fully 
turbulent. 

Variation with altitude of the computed heat transfer coefficient, Ch, for hemispheres 
of various sizes at a flight velocity of 20 km/s is shown in Fig. 3a, and the variation 
for a 30 m (diameter) hemisphere at three different flight velocities is shown in Fig. 
3b. The maximum value of Ch occurs between 55 and 35 km altitude (decreasing 
altitude with decreasing size), and the values appear to collapse around 15 km 
altitude. It should be borne in mind that these results are from computations that 
assume no ablation, and magnitudes will change once material response is included. 
However, preliminary values of heat transfer coefficient are now available as a 


function of velocity, altitude (proxy for density), and size, and the tables can be 
incorporated into the trajectory simulation tool; an appropriate interpolation scheme 
will be required for use in the mass loss equation. 
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(b) Variation with velocity for a 30 m (diameter) hemisphere 

Figure 3. Computed heat transfer coefficient variation with altitude for: (a) hemispheres of 
various diameters for a velocity of 20 km/s, and (b) a 30 m diameter hemisphere at 3 velocities. 


As a test of how far we could push our grid generation and simulation capability and 
get an idea of the level of effort necessary to execute the computations, we have 
also made an attempt at computing the flow around a real shape - asteroid Itokawa. 
We took the available CAD file and scaled the geometry from its original size of »500 


m down to 20 m. This one-off computation was performed for a velocity of 20 km/s 
and a stagnation pressure of 30 bar. The choice of orientation was arbitrary, but 
additional computations can be performed for reasonably small variations in total 
angle of attack. Contours of surface pressure and convective heat flux are shown in 
Fig. 4. Of note in the results is the increased pressure in the neck region of 
geometry - the pressure is roughly 1.5 times the stagnation pressure. We speculate 
that in the case of multi-lobed meteoroids, fragmentation is likely to start in this 
region. 
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(a) Pressure (b) Heat flux 

Figure 4. Surface pressure (in bar) and convective heat flux (in W/cm^) distributions for a 1/38- 
scale model of Asteroid Itokawa. Results are shown for a velocity of 20 km/s and a stagnation 
pressure of 30 bar. The flow is into the plane of the paper. Again, radiative heat flux is not 
included. 


5.3. Additional computations 

The influence of coupling of radiation to the flow was tested for a small number of 
cases (1 m diameter hemisphere at 12 and 16 km/s for a stagnation pressure of 10 
bar). The change in the heat transfer coefficient with coupling was found to be 
between 10-20% lower than without coupling of flow and radiation. Clearly, the 
influence of radiation on the flow and vice versa is significant. Since our coupled 
computations are very slow and resource intensive, alternate methods [16-18] have 
to be explored. 

A good portion of the current effort has been focused on flow computations because 
environments are needed for thermal and thermal-structural responses of materials, 
either as initial conditions, or as boundary conditions. Furthermore, these response 
models require more than just properties measured at room temperature; they 
require the temperature dependence as well. For meteoritic materials of interest to 
the PD effort at NASA ARC, not all properties are available, and these properties 
(and their measurement) are the purview of a PFIA characterization team. In the 
interim, elemental composition of an ordinary chondrite was provided to the entry 
physics team [2]. This composition was used within a suite of materials tools to 
construct what are known as B-prime tables [19]. The B-prime parameter is a non- 
dimensional mass blowing rate typically tabulated over a range of temperatures and 


pressures. B-prime tables are required for the analysis of ablative materials used in 
thermal protection systems of entry capsules. Tables of B-prime for an ordinary 
chondrite over a temperature range of 500-6000 K, and pressures ranging from 0.1 
atm to 1000 atm, were constructed and a computation was performed with Fiat [12] 
at the stagnation point of a test article (made of ordinary chondrite material) at arc jet 
conditions. While these conditions do not correspond to flight, the computations do 
provide a good first test of the model implementation in Fiat. 

Two- and three-dimensional discretizations of a 10 m diameter hemispherical 
asteroid were constructed for the first set of finite element thermal-structural 
analyses in Marc [14]. In these early computations, we have assumed the material 
to be uniformly dense fused silica quartz. For 2D/axisymmetric analysis, 4-noded 
quadratic elements were used, with mesh refinement at the wetted boundary to 
accommodate the expected steep temperature gradient due to high convective and 
radiative heating. The pressure, shear and heat flux from flow computations were 
imposed as boundary conditions on the wetted surface. The back surface was 
constrained to avoid any rigid body motion during the analyses. For 3D analysis, 8- 
noded hexahedral elements were used instead of tetrahedral elements; 8-noded 
elements provide better convergence in Marc. The coupled thermal-structural 
computations for the 2D/axisymmetric geometry showed that the contribution to 
stresses in the body due to thermal gradients and shear were insignificant compared 
with stresses due to pressure loads. Therefore, only structural analysis was 
performed for the 3D geometry subject to just pressure loads. Both the thermal- 
structural and structural computations helped establish grid and resource 
requirements for such analyses. In the longer term, we plan to perform thermal- 
structural and structural analyses using properties that are more appropriate for 
asteroid material, and include cracks and voids in the structure. 

Assuming that a meteoroid is porous, we investigated a hypothesis that the 
penetration of hot plasma from the shocked gas into the porous medium of the 
meteor could cause it to fragment. Using the data available for meteor permeability 
[20], we found that the depth of plasma penetration compared to the size of objects 
we are interested in (10-30 m) is very small. Furthermore, including the effect of 
surface heating impedes the flow of gas into the medium - the gas is cooled as it 
penetrates deeper into the medium. It should be noted that these simulations did not 
account for ablation. For a recessing surface there would be further reduction in the 
depth of gas penetration. The problem needs to investigated further, i.e., with the 
inclusion of surface recession, to draw definite conclusions. 

5.4. Multibody computations 

Our early efforts have been largely focused on computations for a single body, i.e., 
the entry object is a geometrically regular body of uniform mass distribution. In 
reality, the entry object could be a loosely bound (via some cohesive force) collection 
of smaller objects, or an irregular shaped object with two or more lobes. 
Consideration of multi-lobed shapes (regular or irregular) and collection of multiple 
objects is a necessary first step towards understanding not only fragmentation but 
also how multiple bodies interact with each other. For instance, multiple bodies in 
close proximity will experience increased aerodynamic loads and increased local 
heating (convective) due to shock-shock interactions, which in turn could affect the 


light production from the gas cap and wake. For multiple bodies that are sufficiently 
separated in space one would have independent motion of several bodies, which 
could lead to a reduction in light production because of reduced radiating volumes of 
the gas cap and wake. Results from our first attempt at handling lobed shapes and 
multiple bodies are shown in Fig. 5. Since the geometries considered possess a 
rotational axis of symmetry, temperature contours are shown only in the pitch plane 
of symmetry with flow going from left to right. Clearly the shape, size, and number of 
bodies alter the flow structure in the wake. These very recent results have not yet 
been processed with the radiation computation tool Neqair. It should be noted here 
that our current methodology is only capable of handling “static” arrangements of 
multiple bodies, i.e., the flow computations do not include the dynamics of the bodies 
based on the aerodynamic forces (and moments). However, the relative 
arrangements of bodies might tell us something about the influence of shock-shock 
interaction on the light production. 
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Figure 5. Pitch plane contours of temperature (in K) for a freestream condition of 20 km/s and 
stagnation pressure of 30 bar. 



6. Concluding Remarks and Forward Work 

Computations have been performed for hemispherical shape ranging from 1 m 
diameter to 100 m diameter and velocities ranging from 12 km/s to 20 km/s. For 
each combination of diameter and velocity, the altitude has been varied such that 
stagnation pressures between 0.1 bar and 300 bar would be achieved. From the 
computed solutions the heat transfer coefficients were determined for use in a 
trajectory code. Furthermore, we have made initial explorations in: (a) thermal 
response of meteoritic materials, (b) structural and thermal structural response of 
hemispheres subject to entry loads, (c) flow around collections of bodies, and (d) 
flow past a large size irregular object. More details of the present work can be found 
in NASA Technical Memorandum [21]. The present short effort has given the entry 
physics team a good idea of the scale and complexity of the problem of atmospheric 
entries of large size PHAs. However, there is a lot more to do to meet the 
overarching objectives: (i) modeling of energy deposition, and (ii) modeling of 
fragmentation. 


Before taking on the long term objectives the team is currently focused on: 

1. Enhancements of the thermodynamic and transport models to include multiply- 

ionized species - and Expanding the species set will 

open up the velocity space for analysis, so we can consider velocities greater 
than 20 km/s. Relevant thermodynamic properties are available in the work of 
Capitelli et al. [22], Their data have a temperature limit of 50,000 K, and they 
have considered three levels - 250, 500, and 1000 cm'^ - of lowering of 
ionization potential.. 

2. Development of an efficient process to construct synthetic light curves from high- 
fidelity solutions for a single body at multiple points along its flight trajectory. An 
early version of such a process is already in place [23], and needs to be updated 
with the current versions of the simulation tools and utilities. The assumption here 
is that the trajectory simulation tool has been updated with time-varying heat 
transfer coefficients. Even if the high-fidelity solutions do not account for material 
thermal response, the synthetic light curves when compared against actual flight 
observations will provide a good idea of enhancement/reduction in light 
production. The main complication here is fragmentation, which will require a lot 
more work. 

In the longer term, apart from flow simulations, the two main focus areas are: 

1. Material thermal response and its coupling to flow computations. We recognize 
that the problem is complicated by the fact that the flow simulations have to be 
able to handle rapidly moving (recessing) surfaces, and intense radiative heating. 

2. Structural response, primarily from a brittle fracture point of view. Initial 
simulations of thermal-structural response have shown that heat penetration into 
the interior is slow compared to the time of flight, and that only structural 
response needs to be studied. Although this reduces solution turnaround time, 
there is still the issue of voids and cracks in the structure. Furthermore, the 
influence of these on the aerothermal environment is not known. 

Given our current understanding of the processes at work, and uncertainties about 
the material characteristics, including composition and state, our goal is to perform 
simulations starting with simplified models for well-known bolides, and to improve the 
fidelity of models in a systematic way, e.g., inclusion of physical features such as 
cracks, voids, etc. The light curves of observed bolides serve as useful benchmarks 
for comparing results from independent modeling and simulation. However, it is not 
yet clear if light curves tell us something about fragmentation processes at work. In 
some cases infrasonic signatures could provide a clue about major fragmentation 
events. Furthermore, the size distributions and impact footprints of meteorite falls 
could provide valuable guides to fragmentation altitudes. 

Finally, we note that data from bolide observations alone may or may not be 
sufficient. As in the case of building spacecraft entry systems, well designed ground- 
based tests may help in verifying hypotheses we develop during the course of the 
present work, and perhaps help quantify the uncertainties in the predictions. We 
could use arc jets, ballistic range, shock tubes, and laser heating facilities. Even with 
such testing, the issue of traceability of ground-test data to flight will remain, since 
almost all of NASA’s test facilities were intended to address atmospheric entries of 
capsules (or winged vehicles). Clearly, bringing higher fidelity modeling to bear on 


the PHA problem is fraught with challenges, but it is an overdue endeavor that is as 
vital as any other in NASA’s mission to benefit all humankind. 
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